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Dynamics of interacting cold atomic gases have recently become a focus of both experimental 
and theoretical studies. Often cold atom systems show hydrodynamic behavior and support the 
propagation of nonlinear dispersive waves. Although this propagation depends on many details 
of the system, a great insight can be obtained in the rather universal limit of weak nonlinearity, 
dispersion and dissipation (WNDD). In this limit, using a reductive perturbation method we map 
some of the hydrodynamic models relevant to cold atoms to well known chiral one-dimensional 
equations such as KdV, Burgers, KdV-Burgers, and Benjamin-Ono equations. These equations 
have been thoroughly studied in literature. The mapping gives us a simple way to make estimates 
for original hydrodynamic equations and to study the interplay between nonlinearity, dissipation 
and dispersion which are the hallmarks of nonlinear hydrodynamics. 
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I. INTRODUCTION 

Recently, there has been a flurry of experiments that 
study the collective behavior of systems of particles with 
both bosonicPS] anc j fermionicPK^ statistics. There have 
also appeared experimental realizations of cold atom sys- 
tems with long range interatomic interactions (such as 
dipolar boson d 13 " 15 !, dipolar fermions^, ions, and Ryd- 
berg atoms) . Investigation of dynamics of such systems 
is interesting and important for many reasons. The dy- 
namics can provide an important information about the 
nature and strength of particle interactions. In many 
cases the nonlinear dynamical evolution leads to a for- 
mation of shock waves and solitons. Out of equilibrium 
dynamics in the presence of defects and disorder (see, 
e.g., Ref. IT7| can shine light on dissipation and localiza- 
tion phenomena. One of the most intriguing and impor- 
tant direction is the search for universality in dynamical 
properties of systems governed by different microscopic 
Hamiltonians 

The dynamics of cold Bose gases is being extensively 
investigated both theoreticall y ^ 19 * 20 * and experimentally. 
Experiments on Bose-Einstein condensates (BEC) in- 
clude sound velocity measurements^ and collisions of 
atomic clouds which lead to th e obs ervation of dispersive 
shock waves and soliton trains^U. A dissipative trans- 
port in BEC has been recently investigated in Ref. [T7] 
The propagation of matter-wave soliton trains in BECs 
has been realized in Ref. 22 . A strongly interacting limit 
of interacting Id bosons a.k.a. the Tonk's gas analogous 
to free fermions has been realized in Ref. 1231 

The dynamic s of Ferm i systems^ has also been of great 
interest recentlySSEHl. For instance, in Ref. [55] the spin 
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transport was studied by colliding two oppositely spin- 
polarized clouds of fermions. In Ref s . ITU1 and ITT1 some as- 
pects of nonlinear hydrodynamics (such as shock waves) 
in unitary Fermi gases were investigated. These experi- 
ments have also motivated new numerical studies of non- 
linear dynamics in Fermi systems^. 

It is also worth mentioning that in addition to con- 
ventional Fermi and Bose many body systems there ex- 
ists a remarkable family of models interpolating between 
fermionic and bosonic systems. This is the family of 
Calogero models. The exact nonlinear collective descrip- 
tion for these models is knowcPsHUl and has a form of 
equations of hydrodynamic type. These hydrodynamic 
equations are integrable and exhibit features such as soli- 
tons and dispersive shock waves. The equations are 
related to a known integro-differential equation - the 
Benjamin-Ono equatio n 29 " 31 !. Unlike other integrable 
models, the Calogero model and its hydrodynamic de- 
scription retain integrability even in the presence of an 
external harmonic potential. In particular the solutions 
of multi-soliton type have been found for these models in 
Ref. [3J 

Because of the variety of systems and models used 
in studying nonlinear dynamics of cold atom systems 
the understanding of major effects resulting from the in- 
terplay between nonlinearity, dispersion and dissipation 
within a simple unifying description would be very use- 
ful. Fortunately, this description is well known in the 
limit of weak nonlinearity, dispersion, and dissipation. It 
is summarized by the Korteweg-de Vries-Burgers equa- 
tion (KdVB) (see eq. [25] below) . The main goal of this 
paper is to connect this universal picture to some partic- 
ular models used in collective descriptions of cold atom 
systems. For simplicity we concentrate here on simple 
Galilean invariant fluids although many features of non- 
linear dynamics described in this paper can be found in 
more complex systems as well. 

The following is a brief outline of this paper. We start 
by constructing a rather general one-dimensional model 
of a simple Galilean invariant fluid in Sec. [TTJ We lin- 
earize this model in Sec. |III| and proceed to the deriva- 
tion of the effective KdVB using the reductive pertur- 
bation method in Sec. IIV Al We describe the effects of 
nonlinearity, dispersion and dissipation for the effective 
KdVB equation and introduce corresponding scales and 
limits. These scales and limits are presented in the tri- 
angle phase diagram in Fig. [T] We conclude the qualita- 
tive description of KdVB dynamics in Sec. |V| describing 
two different scenaria for shock wave formation in the 
presence of both dissipation and dispersion. Finally, in 
Sec. | VI | we give relations between the coefficients of the 
effective KdVB equation and parameters of some models 
used in cold atom studies. We conclude with some open 
questions and possible generalizations of our results in 
Sec. IVIII and describe non-KdVB universal behaviors for 
some systems with long range interactions in Appendix 

m 
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Figure 1. A schematic phase-diagram summarizing various 
hydrodynamic regimes of KdVB equation. We consider a 
typical step profile shown in the left insert characterized by a 
typical width W and an amplitude U. The directions towards 
the corners of the triangle correspond to the increase of cor- 
responding v, a and ("-terms of the KdVB equation (251 eval- 



uated on the typical profile. Corresponding to the strength of 
those terms the diagram is divided into regions corresponding 
to diffusive, dispersive and shock wave regimes with domi- 
nating dissipation, dispersion and nonlinearity terms, respec- 
tively. The edges of the triangle correspond to the exactly 
solvable equations with one of the terms vanishing. See sec- 
tions [IV]V] for details. 



II. HYDRODYNAMIC MODEL 

The goal of this section is to construct a hydrodynamic 
model (hydrodynamic equations) describing a rather gen- 
eral one-dimensional motion of a simple fluid. An effec- 
tive one-dimensional motion usually describes a class of 
solutions of three-dimensional hydrodynamics obtained 
by a requirement on hydrodynamic fields p{x, y, z; t) — 
p(x, t) etc. or through more complicated procedures of 
a dimensional reduction, e.g., due to the presence of a 
quasi-one-dimcnsional trap in cold atom systems (see 
e.g., Sec. IVlland Refs. [10] and HU . An effective "truly" 



Id model can also appear as a r esult of the transverse 
quantization in a quantum systeirP^LaS] j n ^ ne f n ow i n g 
we assume that the reduction to one dimension is al- 
ready performed and outline other assumptions that we 
use to construct an effective one-dimensional hydrody- 
namic model. 

a. A fluid is one- component. We do not consider 
here more complicated theories of mixtures of fluids. In 
particular, we do not consider the two-fluid hydrodynam- 
ics of superfluids with non- vanishing normal component. 
If a system under consideration is in the superfluid regime 
we assume here that the normal component can be "inte- 
grated out" resulting in a hydrodynamics of an effective 
single-component fluid with some additional terms (such 
as viscous terms) generated by the normal component. 
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b. Entropy generation is very small. This assump- 
tion allows us to consider closed equations for density 
and velocity fields only and assume that the fluid motion 
is isentropic. Combined with the first assumption this 
means that the only relevant degrees of freedom can be 
described by one-dimensional velocity v{x,t) and density 
p{x, t) fields. 

c. Locality. We assume that the energy functional is 
local in density and velocity fields and equations of mo- 
tion contain only local values of fields and their deriva- 
tives. This constraint can be relaxed (see examples in 
Appendix |A| . 

d. Galilean invariance. In the presence of Galilean 
invariance the energy density of the fluid is given by 
pv 2 /2 with other terms either independent of velocity 
or proportional to velocity gradients. Systems without 
this invariance cannot be described by a simple hydrody- 
namic Hamiltonian (JT|). We focus on Galilean invariant 
systems in this work for simplicity but do expect that the 
chiral differential equations of the form ( 25 1 still provide 



a rather universal effective description of collective dy- 
namics even in the absence of Galilean invariance (e.g., 
for the effective description of XXZ spin chain in external 
magnetic field^S). 

e. Energy functional contains density derivatives up 
to the second order. This approximation is known as 
the Boussinesque approximation. In the absence of dissi- 
pation this assumption combined with previous assump- 
tions leads to the following hydrodynamic Hamiltonian 



H = dx 



pv 
~2~ 



pe(p) + A(p) 



{d xP f 
Ap 



(1) 



Here A(p) is some function of density. We notice 
here that if A = const the gradient term becomes the 
Madelung term well-known in superfluid hydrodynam- 
ics. The hydrodynamic Hamiltonian equipped with 
Poisson's brackets 

{p(x),v(x')}=d x 5(x-x') (2) 

generates hydrodynamic equations of an ideal fluicP^Il. 

/. Dissipative function is quadratic in velocity gradi- 
ents. The positivity of the dissipative function does not 
allow for terms linear in velocity gradients and we have 
generically for (Rayleigh) dissipative function 



F=\\ dxn B (p)(d x vf 



(3) 



where t]b(p) is a bulk viscosity of a ID fluid. In the 
presence of dissipation the hydrodynamic system is not 
Hamiltonian and Hamilton equations following from ( 1|2 ) 
should be supplemented by dissipative terms encoded in 
pi . We arrive at the following equations of motion 



dtp 
d t v 



dx(pv) = 0, 
-d x {i] B dxv) 



(d p A)(d x ^pf - A 



Vp 



(5) 



where w = d p (pe(p)) is the specific enthalpy (same as the 
chemical potential at zero temperature) of the fluid and 
the r.h.s. of the second equation is obtained as —-4^. 
The rate of the dissipation of energy is then given by 2F. 

Let us rewrite these equation in terms of density and 
current (momentum density) j = pv. We have 



dtp + d x j = 0, 
d t j + d :c (T + T') = 0, 
j = pv, 



T = pv 
T' = ~rj B d x v 



P - l - P VAd a 



Ad x \np), 



(6) 
(7) 
(8) 

(9) 
(10) 



Here the pressure P(p) is given by d p P — pd p w. The last 
term of ^ is called the "quantum pressure term" in the 
context of quantum fluids. The viscous part T" of the 
stress tensor is linear in the gradient of velocity and is 
obtained from d x T' = SF/Sv^ 

Many models of recent interest can be casted in the 
form of the equations (4j5l. The expressions for w{p), 
A(p) and t]b(p) for some of those models are listed in the 
Table D 

Equations ( 4|5 1 are nonlinear. The nonlinearity results 
in the steepening of density and velocity profiles during 
the evolution and in shock waves. The quantum pressure 
terms (terms containing A(p)) depend on density gradi- 
ents. They typically give rise to dispersive oscillations 
of densities and velocities. The viscosity term (the r.h.s. 
of ^ contains the gradient of velocity field. It results 
in damping. All these terms play rather different roles 
in fluid dynamics and understanding their interplay is of 
great interest both theoretically and experimentally. 



III. LINEARIZATION 



Let us start by studying the system ( i|5 1 in the limit 
of small deviations from the uniform solution p(x) = po. 
v(x) = 0. We consider this solution as a background 
configuration of fields and linearize equations in Sp(x) = 
p{x) — po and Sv(x) — v(x). We obtain 



Sp + pqSv' = 



Po 



%Po Po 



(11) 
(12) 



Here subscript "0" means that the quantity is calculated 
at p = p . Linearized equations give the dispersion equa- 



tion (for 5p ' 



-iut-\-ikx 



etc.) 



1 



= c A k z + -A)fc 4 - 2iv G ujk 2 



(4) where c 2 = d p P\ Po is a linear sound velocity and v$ 



(13) 



>)B() 

2po 



is propo rtio nal to the kinematic bulk viscosity at p = pq. 
Solving (13 1 and expanding up to fc 3 we arrive to 



w(k) « ±ck 1 



A 



4c 2 



H k 2 



iv^k 2 



(14) 
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For the wave propagating to the right we have 



u)(k) w ck — iv^k 2 



Ac 



-1/' 



(15) 



Here we kept only leading terms representing different 
effects (dissipation v ) and dispersion (A ). We consider 
the limit of small dissipation and dispersion and, there- 
fore, neglected the correction to the dispersion quadratic 
in VQ assuming it smaller than the corresponding Aq con- 
tribution. We keep the term A$k 3 which can be compa- 
rable in magnitude to the term iv$k 2 at small (but finite) 
k if the viscosity v§ is very small. The dispersion ( p~5[ ) can 
be reproduced by the following linear equation 



Ut 



cu x 



au x 



is u x 



(16) 



where a = Aq/(Ac). Here one can substitute either Sp or 
Sv instead of u as Sp and Sv are related linearly by (111. 



To obtain an effective equation for the left moving linear 
wave one should just change x — > —x in ( 16 1. 



In the linear approximation an arbitrary initial profile 
of density and velocity is split into right and left moving 
waves. These waves move with velocities ±c respectively 
and slowly disperse and decay. If the amplitude of the 
initial profile is small but finite one should add nonlin- 



ear corrections to chiral (right and left) equations (16 I. 
In addition to this, the nonlinear terms also couple the 
equations for left and right moving waves. However, as 
the right and left profiles pass each other with finite ve- 
locity (with the relative velocity 2c) they interact only 
for a short time. In the limit of weak nonlinearity the 
coupling between equations is small and will not signifi- 
cantly change solutions. Therefore, it can be neg lectecP 
and most important nonlinear corrections should enter 
the chiral wave equation (16 1 itself. 



The goal of the next section is to justify the above ar- 
guments and to derive ( |16[ ) together with the correspond- 
ing nonlinear corrections using the well-known reductive 
perturbation methooP^H^. 



IV. 



KDVB EQUATION VIA REDUCTIVE 
PERTURBATION METHOD 



The linear equation ( 16 ) suggests that in the long wave 



approximation k — > the dissipative v$ term always wins 
over the dispersive a term. This is true unless the vis- 
cosity coefficient vq is small (or zero) for some reason 
and we are interested in the regime of small but finite k. 
Then the comparison of the first and third terms on the 
l.h.s. of (16 1 suggests that we have to treat all deriva- 



tives as small with their relative values determined by the 
scaling dt = 0(e 3 ) and d x — 0(e) where e is some aux- 
iliary "counting" parameter treated as small parameter. 
This scaling was introduced in Ref. 3H m the derivation 
of the Korteweg-de Vries equation for water waves. The 
only modification we make here is that we consider the 
viscosity coefficient v = 0(e) which would allow us to 
have the r.h.s. of (16 1 of the same order as u t and u xxx . 



The aim of this section is to include a weak nonlinearity 
into this scaling. This goal is achieved by using the so- 
called reductive perturbation method (or a proper power 
counting scheme in modern language). 



A. Korteweg de Vries Burgers equation (KdVB) 



Let us introduce the following scaling scheme 

e 2 /(e£,e 3 i), (17) 

e 2 g(e£,e 3 t), (18) 

eVB, (19) 



Sv - 



where 



VB 



Ct — Xq- 



(20) 



The scaling laws e£, e 3 t and er\B have been already moti- 
vated on the basis of linear equation ( 16 ) in the beginning 
of this section. The e 2 scaling of the amplitudes can be 
obtained by comparing the relative values of, say, the 
term u t ~ e 3 u of ( 16 ) with a foreshadowed nonlinear cor- 
rection uu x ~ eu 2 . Then according to a general reductive 



perturbation method the scaling scheme ( 17|19 ) should 
be supplemented by an appropriate perturbation scheme 



/ 
.9 



/ (0) +e 2 / (1) 



(21) 
(22) 



Subst ituting p = p + e 2 /(e£, e 3 t) a nd v = e 2 g(e^, e 3 t) 
with ( 21|22|19 ) into the system ( 4|5 1 we obtain a system 
of coupled equations which can be analyzed in each order 
of e separately. The first non- vanishing order is 0(e 3 ) 
gives the following two equations 



? (o) 
c 



/CO) 



(0), 



eg 



Po 

(0) 



(23) 



A compatibility of these equations determines sound ve- 
locity c giving the well-known thermodynamic relation 



p w Q . 



(24) 



The equation (23 1 gives a linear relation between /w 
and g(°\ 

The next non- vanishing order is 0(e 5 ). It gives 



a t fW+d € (fWgW) = dt(cfV-p gM), 



d t g {0) + d e 



7 (o)' 



2 



Po 



2 Po 



ah 



We use the relation (23 1 to exclude g^ from these equa- 



tions. Then the difference of these two equations gives a 
closed equation for u = f^: 

u t + (uii£ — ctu^e — vu iS, (KdVB) (25) 
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with 



l/Sc 
2 V2p 



WqPO 

2c 



A) 



9c 
9p 



Ao 

4c' 

i/ = 



2ft) ' 



(26) 

(27) 
(28) 



A similar equation for can be obtained from (25 I us- 
ing ( 23 ) . The equation ( 25 ) is known as the Korteweg-de 



Vries-Burgers equation (KdVB). It describes the correct 
far-field behavior of solutions of ( |4|5[ ) in the limit of weak 
nonlincarity, dispersion and dissipation. A linearization 
of (25 1 gives |T6|) as it was expected. 



We note here that in the conventional form of KdVB 
equation (and KdV equation, see below) the a term is 
usually written with plus sign. There is, however, a sim- 
ple transformation 



u{x) 



-x) 



(29) 



which maps the solutions of equations with opposite signs 
of a to each other. Using this transformation we can 
borrow any solutions of conventional KdVB and KdV 



equations to find the corresponding solutions of ( 25 ) 



B. Scales and phase diagram of KdVB 

KdVB has been extensively studied (see e.g., Refs. HSJ- 
|4"9"1 It is a non-integrable equation and only few travel- 
ing wave solutions of KdVB are known analyticalljBH. 
KdVB has a reach dynamics which includes an intricate 
interplay of nonlinearity, dispersion and dissipation. An 
example of numerical analysis of KdVB can be found in 
Ref . HU , . . , 

The effective KdVB equation p5j with fl26|27|28| de- 
scribes the universal limit of weak nonlinearity, disper- 
sion and dissipation of one-dimensional solutions for a 
large class of hydrodynamic systems. Depending on ini- 
tial conditions different terms of ( 25 1 play different role 
at different times of evolution. The goal of this section is 
to summarize different regimes of hydrodynamic behav- 
ior and to show how to make simple estimates for onsets 
of these regimes. 

To understand qualitatively the role of different effects 
let us consider an initial profile u(x, t — 0) which is char- 
acterized just by two scales: a typical width of the profile 
W and a typical amplitude of the initial profile U. One 
can think of this initial profile as of the step-like function 
w(£) shown in the insert of Fig. [T|or as of the single lump 
of density (velocity). We use the parameters U and W 
together with the values of coefficients (, a, v to form 
characteristic times corresponding to the different terms 
of KdVB. We characterize the strength of nonlinearity, 
dispersion and dissipation by the following inverse times 
(characteristic frequencies) 



(uw-\ n c 



aW~ 



vW- 



(30) 



respectively. The relative values of these frequencies de- 
fine the regime of the evolution of the initial profile de- 
scribed by KdVB. Different regimes are summarized in 
the Figure [T] It is also convenient to introduce the spatial 
scales 



v 

Ct7' 



(31) 



These scales are defined so that, e.g., when W ~ W^ v the 
£- and z/-terms of KdVB are of the same order fi^ ~ Q, v . 

The center of the triangle in Figure [l] represents an ini- 
tial profile for which all three scales arc of the same order, 
i.e. nonlinearity, dispersion and dissipation are equally 
important in the beginning of the evolution. Generally, 
different initial profiles correspond to different points of 
the triangle so that the corresponding parameters grow 
from zero at the respective side to the infinity in the re- 
spective vertex of the triangle. 

In the absence of nonlinearity (the bottom part of the 
triangle, VLq <C Q a ,v) the dynamics is approximately lin- 
ear and is described by |T6| ). It is either diffusive (left 
part of the triangle, fl a <C fi v ) or initially dispersive 
(right part, 57„ <C £l a ). The dispersive evolution even 
if dominant in the beginning of the evolution tt a 3> 0„ 
leads to the growth of the width of the initial profile and 
the decrease of gradient terms. The dispersion term will 
become of the same order as the diffusion one at the time 



tn 



■w. 



(32) 



At that time the width of the profile becomes of the or- 
der of W ua (31) and after that moment the dispersion 



becomes subdominant to the diffusion. 

The evolution is much more interesting and compli- 
cated when the nonlinear term dominates for an initial 
profile. In this case the evolution is initially described by 
the Riemann-Hopf (or inviscid Burgers) equation 



u t + (uii£ = 0, (Riemann-Hopf) 



(33) 



which can be easily solved for any initial profile u(x, t = 
0) = f(x) giving u(x,t) implicitly, as a solution of u = 
f(x — Qut), where the unknown u enters both left and 
right hand sides. For our typical initial [/TV-profile this 
solution is well defined at small times becoming multiply- 
valued for times after 



tr 



n- 1 = w/(cu). 



(34) 



The time t c is known as the time of "gradient catastro- 
phe" and is defined as the time at which u x becomes 
infinite. The classical problem (33) is ill-posed at larger 



times and has to be regularized by higher gradient cor- 
rections which in our case are either dispersive (a-term 
of ( 25 1 ) or diffusive (u-term). According to the relative 
strength of subdominant diffusive and dispersive terms 
in the regime of strong nonlinearity (the upper part of 
the triangle phase diagram) one has either formation of 
classical dissipative shock waves (left part) or dispersive 
shock waves (right part). We discuss these shock wave 
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regimes separately in SecfV] Before going to this discus- 



sion we consider important limits of KdVB equation ( 25 ) 



corresponding to the left and right sides of the triangle 
phase diagram of Fig. [T] 



C. Dissipative limit: Burgers equation 

If the dispersion does not play a significant role in the 
evolution (Tl a <C Qv) one can neglect a-term in (25 1 and 
arrive to the well-known Burgers equation 



u t + (uu£ — uu^c (Burgers). 



(35) 



The regime described by ( 35 1 corresponds to the left side 
of the triangle phase diagram of Fig. [T] The different 
points of the left side the triangle correspond to different 
relative values of nonlinear and dissipative terms of ( 35 1 . 



We note here that one can arrive to the Burgers equa- 
tion (35 1 more formally starting directly from (4|5 1. In- 
stead of (19 1 one should take t\b — 0(e°). Then the 
scaling scheme ( 17|18 1 is not consistent and should be 



replaced by the Burgers scheme^ 

Sp = ef{ex,e 2 t), 
v = eg(ex, e 2 t) 

with the corresponding perturbative expansion 



/ = / (0) 
9 = 9 i0) 



(36) 
(37) 



(38) 
(39) 



instead of ( 21|22[ ). Expanding in e we arrive to the re 



lation ( 23 1 and to the Burgers equation (351 as a close d 
equation for u = /(°) with parameters given by ( 26|28 1. 
The Burgers equation ( 35 1 is often used as a 



model describing classical dissipative shock waves which 
result from the interplay between nonlinearity and 
dissipation.^ It can b e analytically solved via Cole-Hopf 
transformation u = — log</> which results in a diffu- 



sion equation (f) t — v(j) x 



50 51 



There is a simple exact solution of ( 35 1 given by a 



step-like profile traveling with a constant velocity V to 
the right 

u{i,t) = -(l-tanh^-^C-W)]). (40) 

Here the parameters of the step are related to the velocity 
V as 



W 

u 



2u/V 
2V/C- 



AW, 



(41) 
(42) 



The width of the solution is given by (41 1, i.e., essentially 
by the scale W^ v . In the limit v — > this width goes 
to zero and the solution ( |40| d escribes a discontinuous 
shock front. The solution (|40[) is defined when Am = 

U^-oo - U£-^+oo = U > 0. 



D. The limit of no dissipation: KdV equation 



In the limit when dissipation is absent or very small 



at a given wavelength the v term of KdVB ( 25 1 can be 
dropped and we arrive to the celebrated Korteweg-de 
Vries equation (KdV) 



u t + (uiif: — au^ = (KdV) 



(43) 



with parameters determined by (|26]27j). It is a purely 
dispersive equation. It is integrabl d ^ * 53 ! and possesses 
infinitely many conserved quantities. The first three of 
these integrals are given explicitly by 



h 
h 



2 ' 



C— + cm. 



(44) 
(45) 
(46) 



The integrals /o,i,2 are related to the total number of 
particles, total momentum and total energy of the system 
( 4p ) as follows 

(47) 
(48) 

(49) 



N ■ 



N 
P 



2c 

Po 

c 



h 



do, 
2c 2 



E-Eo = —h + h + w I - 



Pn 



Po 



The higher integrals of motion of KdV are related to more 
complicated symmetries of the problem and are usually 
destroyed by small corrections to KdV which destroy in- 
tegrability. We do not need their exact form in the fol- 
lowing discussion. The conserved quantities ( 47|48 1 and 
( 49 1 on the other hand are related to fundamental space- 



time-gauge symmetries and, therefore, play an important 
role in the hydrodynamic approach. 

The KdV has a solitary wave solution (soliton) moving 
to the left with velocity V^l 



u(£,t) = -C/cosh" 2 (W^it + Vt) 



(50) 



This solution corresponds to a local depletion of the par- 



ticle density (minus sign in (50)) and is known as the 



dark soliton. The width W of the soliton and the ampli- 
tude U of the depletion in ( 50 1 are defined by its velocity 
V as 



W = {Aa/V) 1/2 = 
U = 3V7C • 



12 W Ca , 



(51) 
(52) 



While the parameters of the soliton ( 50 1 are given ex- 
actly by ( 51|52 l they could also be estimated from the 
condition that the nonlinearity and dispersion scales of 
the soliton solution should be of the same order rt a ~ 
(exactly we have: D,^ = 12f2 Q )J^ 



The solution ( 50 1 corresponds to the special point on 
the right side of the triangle phase diagram in Fig. [T| 
These solutions (as the effective KdV itself) are to be 
trusted only when the relative depletion U /po ~ V/c is 
small, i.e., V«c. 

The total number of depleted particles in the soliton 
solution Q50|51|52 l is given by the zeroth integral of mo- 
tion dUb 



U(x,t) 



u(£,i)d£ 



-2UW 



(53) 



The values of other integrals of motion on the KdV soli- 
ton solution ( 50 1 are 



2 , VVaV 
h = tU 2 W = 12- 



V 

=- ? », 



15 



CU 3 W 



6 V 2 



(54) 
(55) 



C 2 

72 V 2 

\l aV = 

5 C 2 5 C 

Using above expressions and ( |47|48|49| we obtain the 
dispersion of the soliton 

(P — cn) 2 

E(P) -E - nw =c(P-cn)+ y - ; (56) 

2m* 

with the "effective mass" m* given by 

ben 

m = 

The form ( [56] ) is similar to the dispersion of particles with 
quadratic spectrum except that the effective mass ( 57 ) is 
velocity (and momentum) dependent (see eq. 53 ) . 



(57) 



Another important exact solution of KdV is the peri- 
odic traveling wave solution given by 



u(£,t) = -£/cn 2 (w -1 (£ + Vt) rnj 



(58) 



In Eq. [58] cn(y|m) is the Jacobi elliptic function of mod- 
ulus m (0 < m < 1). The modulus defines the period of 
the solution in £ (the period L itself is L = W F(%/2, m)). 
In the limit m — > 1 the period L — > oo and Jacobi ellip- 
tic function reduces to a periodic array of well-separated 
solitons (50 1. 



So far we discussed very special solutions of KdV corre- 
sponding to the "Solitons" point in Fig. [I] What will hap- 
pen to the initial profile with typical dimensions such that 
the nonlinear term is dominant il a ? In this case 

the profile will initially evolve according to the Riemann- 
Hopf equation (33 1. However, at the gradient catastrophe 



time ( 34 1 large gradients will develop and the dispersive 
term will become of the order of the nonlinear term. De- 
tails of the evolution after that point do depend on a 
particular shape of the initial profile. However, generi- 
cally, the modulated periodic solutions of KdV will be 
generated providing oscillating (a.k.a. dispersive) shock 
fronts^. This phenomenon is referred to as Dispersive 
Shock Waves (DSW) in contrast to conventional or dissi- 
pative shock waves which occur in equations of Burgers 
type. Without dissipation the steep profiles described 
by KdV generate oscillations and eventually decay into 
trains of well-separated solitons. 




Figure 2. The evolution of the initial step-like profile gov- 
erned by the KdVB equation ( 25 1 in the dispersive shock wave 



regime. 



SHOCK WAVES IN KDVB 



So far we have discussed the dynamical regimes of 
KdVB corresponding to the sides of a triangle phase di- 
agram Fig. [T] These regimes are realized when one of 
three scales of KdVB ( 30 1 can be neglected compared to 
two others. If all three scales are significant the different 



ordering of these scales roughly divide the phase diagram 
of Fig. [T] into six parts. Up to this point we treated three 
scales as equal but to understand the KdVB dynamics 
of a typical profile even better one should realize that 
the the nonlinearity, dissipation and dispersion enter the 
game in a very asymmetric way. Indeed, we have already 
seen that in the case when the nonlinearity dominates 
an initial evolution (f2^ is the largest scale) the evolution 
itself will result in large gradients in finite time which is 
of the order of the gradient catastrophe time (34 1. At 



that point either dissipation or dispersion will become of 
the order of nonlinearity and our discussion of further 
evolution should start from the different part of the tri- 
angle phase diagrarrPH Another example of this type is 
given by the linearized KdVB (16 1 with dominating dis- 



persive term (f2 a 3> £l v ). The corresponding initial pro- 
file will disperse almost non-dissipatively up to the time 
of the order of t au ( |32| at which point the dissipation 
term will become of the order of the dispersion term and 
the further evolution will be diffusion-like. In contrast 
to these examples, if the dissipation term is dominant in 
the beginning it remains dominant and the evolution is 
diffusion- like at all times. 

In the following we apply the understanding of the 
asymmetry of different KdVB terms to describe two dis- 
tinct shock wave regimes of KdVB. We consider an ini- 
tial profile of the step form shown in the insert of Fig. [T] 
Shock waves appear in the regime when nonlinearity is 
the most important scale in the problem, e.g., if the width 
of the initial profile W is sufficiently bigger than Wq v and 
Wfa plj ). Then the nature of shock waves depends on 
the relative strength of dispersive and dissipative terms. 

The goal of this section is to describe the major fea- 
tures of classical and dispersive shock waves and show 
how to estimate typical scales of these shock wave solu- 
tions. 
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A. Classical shock waves (CSW) 

To have a clear separation of scales we consider the 
following relation between the width of the initial step- 
like profile and scales of KdVB (31) 



w» w Cv » w Ca , w ua . 



(59) 



which means ^ £1„ ^> ^/Tt^Tl^. In this case the dom- 
inant term of KdVB for the initial profile is the nonlinear 
£-term and the evolution of this profile is determined by 
the competition between dominant nonlinear and sub- 
dominant dissipative i/-term. 

Up to the time of the order of t c (34 1 the nonlinear 



term will dominate the evolution which will be described 
by the Riemann-Hopf equation ( 33 1 . Then the decreas- 
ing width of the profile (the width of the wave front) will 
reach the scale W^ v ( 30 I and an approximate balance be- 
tween nonlinearity anadissipation will be achieved. The 
profile will become stationary and will have a width of 
the order of 



ir, 



CSW 



(60) 



The quasi-stationary profile can be roughly described by 
the solution (40 1. The ratio of the dispersion scale to 



the dissipative scale for this profile is given by Q a /Q„ = 
W va /W and will remain very small for all times due to 
(59 1. Therefore, we expect that the dispersive effects are 



not important in the regime ( 59 I and might result only in 



small oscillations on top of the stationary classical shock 
wave solution ( 40 ) . 



B. Dispersive shock waves (DSW) 

The case of the initial step-like profile characterized by 
the width 



w, w va » W Ca » W Cl , 



(61) 



is, probably, the most interesting. The inequalities (61) 



mean that the evolution of the initial profile is defined 
in the beginning by an interplay between the dominant 
nonlinear and subdominant dispersion terms with dis- 
sipation playing some role only at very large times. As 
both v- and a-terms are subdominant to the nonlinearity 
the scale W va does not play any role in the beginning of 
the evolution and its relation to the initial profile width 
W is not very important. We note also that the second 
of the inequalities ( [61] ) is a consequence of the first one 
and the smallest scale Wq v does not play a major role in 
this regime. 

As in the case of classical shock waves for t < t c the 



evolution is governed by the Riemann-Hopf equation ( 33 ) 



and the width of the step decreases reaching W^ a in time 
of the order of t c (34 1. At this time oscillations develop 
at the trailing edge of the step profile 58 These oscillations 
grow in amplitude with the largest amplitude becoming 
of the order of the size of the step U . At this point the 



typical wavelength of oscillations is given by the scale 
W^ a from (31 1. The number of oscillations and the spa- 
tial extent ofthe oscillating part of the shock front con- 
tinue to grow. In the absence of dissipation (KdV) this 
growth continues forever with oscillations evolving into 
the train of well-separated solitons 5 -^. 

In the case of non-vanishing dissipation the front of 
DSW becomes stationary at large times. The width of 
the stationary profile can be estimated in the following 
simple way. As the width of the oscillating shock profile is 
determined by the dissipation we evaluate the magnitude 
of the dissipative term f2„ ( 30 1 at the gradients developed 

due to oscillations, 



i.e. 



W f . 



ir 2 



= ^(U. On 



V(c 'Co. 

a stationary profile this scale should be of the order of 
the nonlinear scale evaluated at the overall width of the 

Wnsw ' 



profile Wdsw, i.e. i\ 

Wnsw 

scales we obtain an estimate 



Equating these 



Wdsw 



a 
v 



(62) 



The formation of the stationary shock front profile of 
DSW takes time of the order of tosw ~ W ua / {CJJ) and 
the number of oscillations in the stationary profile can 
be estimated as 

1/2 

\ - ;;r - [ v ] » 1. (63) 



This condition can be considered as a necessary condition 
of observing oscillations in the stationary shock wave pro- 
file. 

We have to remark here that the presented picture cap- 
tures only major scales of the dispersive shock wave for- 
mation. For example, while we assumed for our estimates 
that the typical amplitude of oscillations is U, the oscilla- 
tions at all amplitudes ranging from to ~ U are present 
in the shock front. As a result even after the formation 
of the stationary shock wave profile for major scale U the 
small amplitude oscillations keep propagating at the lead- 
ing edge of the shock wave profile. The reader is referred 
to the seminal Ref. [55] for details. In Ref. I5"9l Gurevich 
and Pitaevskii used the large parameter separating the 
scales WvajW^a 3> 1 to describe the formation of the 
DSW profile of KdVB analytically using the Whitham 
modulation theory. They described the oscillating part 
of the profile by a modulated periodic solution of KdV 
( p58| (see also Ref. UHlfor a recent discussion of DSW in 
KdV and cold atom dynamics). 

We conclude that in striking contrast to the classical 
shock wave the dispersive shock wave profile has an inter- 
nal structure - oscillations with the typical wavelength 
W^ a and that while the width of the conventional shock 
wave is proportional to v (60), the overall width of DSW 
is proportional to v^ 1 (62 1. 



The main steps of the formation of DSW are illustrated 
in Fig. [2] Although the KdVB equation itself is consid- 
ered in this work only as an approximation to more pre- 
cise hydrodynamics ( 4|5 ) it is important that it gives a 
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very good qualitative and often even quantitative under- 
standing of the latter. Below, in section [VI B| using the 
numerical method of smoothed particle hydrodynamics, 
we present the results of numeric solutions of hydrody- 
namic equations written for a cold Unitary Fermi gas (see 
Fig. |3|. The similarity between these solutions and the 
dispersive shock waves of KdVB described in this section 
are evident. Using the results given in this section we will 
make estimates (in Sec VI B I for dispersive shock wave 



scales in a Unitary Fermi gas system for gas parameters 
similar to the recent experiment}^. 



be significant. In simplest cases these effects can be in- 
corporated into the effective one-dimensional hydrody- 
namic description by adding the Burgers terms vu^ with 
v treated MS M phenomenological parameter. 



Weak coupling (high density) limit 7< 1 



It is well-known that a collective description of (64 1 



in the high density limit 7 <C 1 is given by the Gross 
Pitaevski equation^ (GPE), (see, e.g., Refs. 61. and 



VI. INTERACTING COLD ATOMS 



ihd t ip(x,t)= l- — d xx +g\^{x,t)\ 2 \^{x,t). (66) 



In this section we consider several one-dimensional 
models which recently attracted a lot of interest in con- 
nection with cold atom systems. These are the systems of 
bosons with contact interaction in both weak and strong 
coupling limits and the Fermi gas at unitarity. These 
models do respect the Galilean invariance and therefore 
can be described by the general form of ( 4|5 1 in the hy- 
drodynamic regime. The goal of this section is to relate 
the chemical potential w(p) and the values of coefficients 
A, t\b in f5| to the parameters ( 26|27|28 l of the effec- 
tive KdVB ( p5| description. We delegate the analogous 
discussion of several models with long range interactions 
to the Appendix [XJ The results are summarized in the 
Table D 



A. ID bosons with contact interaction 

One-dimensional bosons with contact interaction can 
be described by the Lieb-Liniger model 60 



H 



h 2 A d 2 



E 



2m ^— ' dx 2 



+ g^2S{xi - Xj) . (64) 

i<j 



The background density of bosons p and the coupling 
constant g define the dimensionless coupling 7 as 



7 



m 

h 2 p 



(65) 



The value of 7 determines whether the system is in a 
weak coupling (7 <C 1) or a strong coupling (7^1) 
regime. 

The model (64 1 is integrable by the Bethe Ansatz 
method for all values of the coupling 7. In particular, 
a general form of w(p) can be found implicitly through 
the solution of the Bethe Ansatz integral equations^. 
For a general value of 7 this solution can be found only 
numerically. Here we discuss only the limits of weak and 
strong coupling in which analytical formulas can be ob- 
tained using the expansions in 7 and 1/7, respectively. 

As the model (64 1 is integrable it does not include 
any dissipation mechanism. However, in more realistic 
modelings of experiments the dissipative effects might 



Using a "hydrodynamic" change of variables, 



(67) 



the GPE can be casted in the form ( 4j5 1 with w, A and 
r] B given by 



w{p) 
A 

VB = 



— P, 

m 

2m 2 
0. 



Then from ( |24|26|27[ ) we have 

hp 



C 



m 

h 3 

m 4 
h 



V7. 



(68) 

(69) 
(70) 



(71) 
(72) 
(73) 



As mentioned earlier a small amount of dissipation that 
arises experimentally can be taken into account by intro- 
ducing v as a phenomenological parameter. Then, the 
dynamics of the model is described by the KdVB equa- 
tion (251 with the values of £ and a. given by ( 72|73 1, 
respectively, and with v as a phenomenological parame- 
ter. In the limit of no dissipation one can describe this 
system by the KdV equation (i.e., KdVB with v = 0). 

The parameters of the soliton solution of an effective 
KdV equation in this case can be found from ( 51 152|53 > 



W = 



U 




and the effective mass (57 1 is 



20 



rn 



(74) 
(75) 
(76) 

(77) 
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Table I. Hydro dynamic descriptions for some cold atom systems. Mapping of the parameters in hydrodynamic equations | 4|5 1 
to the coefficients of the effective KdVB equation. 



Cold Atomic Model 
Chiral differential equation 



w(p) 



Pow 



* _ 3c 

I w 'po 

Z 4c 



Dense Lieb-Liniger Gas 



KdV (Eq. 431 



gp/m 



1/2 



0" 



Pa , 



—^2 2 2 



Tonk's Gas 



Riemann-Hopf (Eq. 33 1 



2m 2 



0" 



Plane- Wave Unitary Gas 



(1 + ^(3^ 



1/2 



7 «'() 



TiWB (Eq. 251 



2 _c_ 

3 Po 



8m 2 c 



Quasi-1D Unitary Gas 
KdVB (Eq. 



25) 



Calogero Gas 
Benjamin-Ono (Eq. 



Quasi- ID Dipolar BEC 



A6I 



l + /5)[f^] ! 9/20 
A 2 /2 



2m 2 ; 2 



2)11 



ml' 2 



1/2 



3c 
r 'P() 



3 _£_ 

4 PO 



HBO 



2m 2 c 



non-local KdV (Eq. A20l 



a In experiments, many atomic gases which to a good approximation behave like Lieb-Liniger gas might still have mechanisms 
of dissipation. The dissipation can be included phenomenologically by adding uu^ to the right hand side of effective 
differential equation which makes it KdVB (dense limit) or Burgers (Tonk's limit) equation. 
6 See the discussion of viscosity in Unitary Fermi gas in Sec. VI B for more details. 



This soliton solution of KdV coincides with the soliton 



solution of GPE equation ( 66 1 in the limit of weak nonlin- 



earity. The latter is known to correspond to a q uasi-hole 
excitation of the quantum Lieb-Liniger model (|64| 63 64 
The soliton solutions should be trusted only in the limit 
|ra| S> 1 and V <C c. The first condition corresponds 
to the applicability of GPE while the second one allows 
us to replace GPE by KdV corresponding to a weakly 
nonlinear limit of GPE. 



2. Strong coupling (low density) limit 7> 1 

In the limit 7 > 1 the che mical po tential is given by 
the following expansion in 1 / ^ 6 ° l 65 l 66 l 



w{p) 



n 2 h 2 p 2 
2m 2 



1- — + ^ + 
-37 7^ 



20 64 V 15 



/y3 



(78) 

where j(p) is understood as a function of the local density 
p (cf. Eq. |65b 



~r \ m 
n p 



(79) 



Using ( 24|26 ) and replacing 7(^0) by its background value 
7 = 7(A)) (65 1 we obtain 



h 

c = — 7rp 
m 

C = — 7T 

m 



1 - 
1 - 



1 

7 
3 

27~ 



(80) 
(81) 



where for the sake of brevity we kept only first couple of 
terms of the corresponding expansions. 

A derivation of the dispersive term A{p) requires some- 
what more involved analysis and will be considered else- 
where. The value A vanishes in the limit 7 —¥ 00 (see 
below) . 

The limit of an extremely strong coupling 7 — > 00 is 
known as the Tonk's limit. This limit corresponds to 
impenetrable bosons which in turn can be mapped to free 
one-dimensional fermions (see Refs. and [37] for recent 
discussions). The collective description in this limit is 
essentially given by A(p) = 0, r\B — and 



7r 2 n 2 p 2 

2m 2 



w(p) 



c = —irpo . 
m 



(82) 
(83) 



The effective KdVB equation becomes the Riemann-Hopf 
equation (Eq. 33 I with 



771 



(84) 



In fact, for ID fermions the left and right moving waves 



are decoupled and the equation (Eq. 33 1 follows directly 



and exactly from the Euler and continuity equations with 
pl. 



An i mportan t point about the Tonks limit lies in the 
absence^ 9 * 68 * 69 ! of the quantum pressure or the zero point 
fluctuation term (i.e., A = 0, a = 0). Therefore, in the 
Tonks limit one does not expect to see the dispersive 
KdV- like oscillations. To account for a dissipative mech- 
anism that might be present in experiments, the dissi- 
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B. Fermi gas at unitarity 




x(iim) 



Figure 3. Num erical dens ity evolution (based on smoothed- 
particle approac h 1 10 * 11 * 70 * ^ ) of a Gaussian initial profile (with 
profile scales similar to the experiment in Ref. 10) in a Unitary 
Fermi gas in cylindrical geometry (see Sec |VI B~3| for details). 
We use the symmetry of the plot and show only its positive 
x-axis part for visual clarity. The red (dashed) plot depicts a 
sit uation when the dissipative term dominates the dispersive 
on e 1 10 * 11 ! (a^ = 1) and the black plot represents a situation 
when both terms are important (a v — 0.1). Albeit in highly 
nonlinear regimes the plots shown can be qualitatively de- 
scribed by KdVB physics of Fig. [T] Indeed, the black plot is 
similar to the one of Fig. [2]while the red (dashed) one is rem- 
iniscent of the evolution governed by the Burgers equation. 
See Sec. IVIBl for related discussions. 



pative term term vu^ can be added with the effective 
equation becoming the Burgers equation ( 35 1 . 



Recently there have been numerous experiments study- 
ing the dynamics of strongly interacting cold Fer mi atoms 
and, in particular, of the Fermi gas at unitarit y! 5 * 26 * 72 * 73 ! 
The Fermi gas at unitarity is highly hydrodynamic even 
deeply in the nonlinear regimeM Interacting fermions 
can also be studied in a quasi-one-dimensional regime 
where the effective ID hydrodynamics is obtained as a 
result of the di mensi onal reduction of the 3D hydrody- 
namic equation a 10 ! 74 *. 

In this section we consider the superfluid hydrodynam- 
ics of the Fermi gas at unitarity. We assume that the tem- 
perature is low so that the gas can be effectively described 
as a simple fluid instead of two-fluid hydrodynamics rel- 
evant for superfluids at finite temperatures. We do take 
the presence of the normal component of the fluid into ac- 
count only through the effective shear viscosity which is 
of the order of r\ s hear ~ ^P with the coefficient of the or- 
der of one under typical experimental conditions^. The 
bulk viscosity should be absent in the hydrodynamics of 
the gas due to the scaling invariance at unitarity^. 

In the following we consider two types of one- 
dimensional behavior of nonlinear waves in unitary Fermi 
gases. The first one is the plane wave propagation 
through the 3D Fermi gas (considered in Sec. VI B 1 1 and 
the other one is the propagation of waves in the Fermi 
gas confined to elongated cigar-shaped traps where the 
effec tive ID reduction of 3D hydrodynamics is usecP^ (see 
Sec. VIB2I. The results are summarized in the Table [T] 



1. Plane waves in 3D unitary gas 

The 3D superfluid hydrodynamics of a unitary Fermi 
gas is given by the continuity and the Euler equations 
with the chemical potential fixed by the scaling invari- 
ance as 

h 2 



w{p) = (l + p) 



2m 2 



(85) 



where j3 w — 0.61 is the Bertsch parameter^SHU] Here the 
density is the 3D density of the Fermi gas. In addition, we 
take into account the gradient corrections coming from 
the Madelung's pressure term withal 



2m 2 



(86) 



and the effective shear viscosity r] s } lear which appear as 
a result of the presence of the normal component. 

We consider the plane wave solutions of the corre- 
sponding hydrodynamic equations assuming that the 
fluid moves along z direction and 



p{x,y,z) 

v z (x,y,z) 



P{z), 
v(z). 



(87) 



The substitution of ( 87 1 into 3D hydrodynamic equations 
results in Id equations ( 4p | with the effective ID bulk 
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viscosity t\b originating from the shear viscosity of the 
3D fluid 



with the transverse oscillator length given by the radial 
trap frequency lj±_ 



°f]B — T^lshear- 



We obtain the effective parameters of KdVB as 



1 + /3 



1/2 



in 



-(3tt po 



and 



c = 2 - 

3 Po 

a= — 
8m 2 c ' 

2 h 

v = o a i" ■ 

3 m 



(89) 

(90) 

(91) 
(92) 



We notice here that the effective dissipation parameter 
v is defined by the shear viscosity of 3D unitary gas 



Vshe 



l n m 



Po- 



Here a v is the dimensionless shear 



viscosity. Putting a n to zero would lead to dispersive 
shock waves (see e.g., Ref. EH). However, the shear vis- 
cosity of the unitary gas has a non-trivial temperature 
behavior! 80 * 81 ^ It has a universal dependence a v ~ T 3 / 2 at 
large temperatures and diverges as a v ~ T -5 at very low 
temperatures. The latter divergence is related to the di- 
vergence of a phonon mean free path (determined by the 
rate of phonon-phonon collisions) at low temperatures.^! 
As a result the shear viscosity of a unitary gas is ex- 
pected to develop a minimum at the temperature close 
to the t emper ature of a superfluid phase transition with 



a 



Substitutin g (90|91|92 l into the criterion for dispersive 
shock waves (63) we obtain a v <C (3C//(16po)) 1 ^ 2 which 
is not possible even extrapolating the amplitude of the 
wave into deeply nonlinear regime U ~ po- We will see 
below that the unitary gas in a quasi- ID trap gives a 
more feasible way to observe dispersive shock waves. 



2. 3d Unitary gas in quasi-ld trap 

There is another important case when the dynam- 
ics of a unitary Fermi gas can be reduced to the one- 
dimensional. This is the case of a unitary gas confined to 
elongated cigar-shaped traps. In this case the longitudi- 
nal waves along the trap can be described by the effective 
reduced ID hydrodynamics^. It is convenient to use the 
effective one-dimensional density p{x, t) equal to the 3D 
density integrated over the transverse section of the trap. 
Then the effective hydrodynamics can be casted in the 
form of Eq. [J] and Eq. [5] where^EH 



w(p) 



2rnH\ 



(1 



/ 157T p l± 



2/5 



(93) 



l± = 



(94) 



The Madelung term of (p| is given by the reduction of 
the bulk Madelung term ( 86 ) and turns out to b<P 



A{ P ) 



9 h 2 

20 77? 



(95) 



The effective one-dimensional bulk viscosity coefficient 
t]b originates from the normal component of the super- 
fluid. It should be obtained in the process of averaging 
of the effective shear viscosity of the 3D superfluid uni- 
tary gas. It is also affected by the boundary of the trap 
where hydrodynamic approximation does not work and a 
more complete kinetic theory should be considered. This 
derivation is beyond the scope of this paper and here 
we treat t\b as a phenomenological parameter. At the 
typical conditions of the existing experiments^ it is 



m 



(96) 



where is a dimensionless parameter of the order of 1 . 
It has been estimated in Ref. [TO] as being in the range 
1 — 10 for the conditions of Ref. 

It was shown in Ref. [TO] that the hydrodynamics ( i|5 1 
with ( |93|96 l gives a very good quantitative description 
of the cloud collision experiment. The dispersive term 
was neglected in Ref. [10] as it turns out to be inefficient 
under the experimental conditions (see below). The hy- 
drodynamic equations can not be solved analytically and 
numerical solutions were used in Ref. 1101 In the following 
we use instead the mapping to the effective KdVB equa- 
tion to identify important scales in the problem and to 
analyze the feasibility of observing the dispersive effects 
in cold Fermi atoms. 

Having the unitary gas dynamics in the form ( 4|5 1 we 
immediately derive the parameters of the effective KdVB 



equation ( 25 I in the limit of weak nonlinearity, dispersion 



and dissipation. The sound velocity is given by (24 1 and 
is equal to 



1+(3\ 1/2 h (\fm poll. x 1/5 



Then ( |26|27|28| give 



mZa V 2 1 + /3 



3c 
9h 2 



80m 2 c ' 
h 



v = ot. 



m 



(97) 

(98) 

(99) 
(100) 
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3. Shock waves in unitary gas in quasi- Id trap 

Let us now use the effective KdVB equation for the uni- 
tary Fermi gas in the elongated harmonic trap to under- 
stand the qualitative behavior and relevant scales for par- 
ticular experimental conditions of Ref. ITUl In the latter 
experiment the two-component Lithium gas was cooled 
in the trap with the radial frequency cjj_ = 2n x 437-ffz. 



The corresponding transverse oscillator length (94 1 can 
be easily found as 



l± = 1.8/xm , 



(101) 



using win = 1.15 x 10 kg corresponding to h/mn ~ 
9.2/im 2 /ms. The one-dimensional density in the middle 
of the trap is 



Po 



1.1 x lO^m -1 , 



(102) 



and could be determined from the equation of the state 
of the unitary Fermi gas, the trap parameters and the to- 
tal number of particles (N ~ 2 x 10 5 ). The temperature 
in the experiment was very low and will be considered as 
zero here (except in using the effective viscosity parame- 
ter). 

The experimental studied the collision of two atomic 
clouds created within the trap using the blue detuned 
laser beam. Here, instead we make our estimates for 
a weakly nonlinear, dispersive and dissipative limit pre- 
sumably described by the effective KdVB equation (25) 



with the speed of sound given by ( 97 1 and other param- 
eters by ( |98|99|100| ) as 



14.5 



/im 

ms 



x 10 



0.65 



x 9.1 



with a r 



10. 



(103) 
(104) 
(105) 
(106) 



Let us now assume that we created a step-lik e pro file 
of the density of the typical scale of U ~ ep or ( 102 1 



U 



e x 1.1 x 10 3 ^to 1 



(107) 



where we keep e< 1 as a parameter. Then the relevant 
spatial scales d3l| are estimated as 



_1 1/zm, 



1/2 



10 /im, 

io~V"i. 



We find that the condition ( 63 1 becomes 



N = 



e^a- 1 » 1 



(108) 
(109) 
(110) 



(111) 



and is definitely not valid for a v of the order of one or 
larger. Therefore, unless the dimensionless viscosity a v 



is somehow experimentally suppressed, the shock waves 
have the conventional dissipative character and disper- 
sive shock waves should not be observed. 

Indeed, the steep density profiles in the experiment of 
Ref. [10] were identified as the dissipative shock waves. 
Stretching our estimates to e ~ 1 to take into account 
the highly nonlinear character of the cloud collision in the 
experiment (of course, KdVB derived in the limit of weak 
nonlincarity can serve here only as a tool for estimates) 
we obtain that N ~ 0.1 and the approximate width of 
the dissipative shock front should be of the order of 



W C sw ~ W c „ ~ lQ/im 



(112) 



This number is indeed of the order of the one observed 
experimentally-^. 

To see whether it is feasible at all to observe the dis- 
persive shock waves in cold Fermi atoms we make the 
following crude estimate. We start with the condition 
( 63 1 rewritten as 



<€.aCU 



and use for the estimate a 
U ~ po- We have 



m 2 c ' 



'I m 1 ^ 



•c l 



(113) 



— and 

Pa 



(114) 



as a necessary condition for the observation of the oscil- 
lations in the stationary shock wave profile. Is it possible 
to lower the effective viscosity to achieve these conditions 
in the experiment similar to Ref. 10? A straightforward 
estimate of the mean free path for phonon-phonon colli- 
sions for conditions of Ref. 10 gives 



A - (T F /T) 9 4 x 10 



-s 



(115) 



We used here the expression^ A 



he 
E F 



x 2.8 x 10~ 5 (1 + 



fif (Jf) with7> w l.l/iKandc= v F /V3 « 19/xm/ms 
calculated at the center of the trap corresponding to a 3D 
density 10 n^o = 6.1/im -3 . 

From (1115 1 we see that the mean free path evaluated 



at the center of the trap becomes of the order of the trap 
size (~ 200 fxm from center to the edge of the trap) al- 
ready at T w O.lTp which is roughly the temperature 
of the cold atom system in Ref. [10] Therefore, we ex- 
pect that if one lowers the temperature the mean free 
path will saturate at the size of the trap and the shear 
viscosity of 3D superfluid will be dominated by the den- 
sity of the normal component of the superfluid p n ~ T 4 
(from kinetic theory rj s hear = PnP^> where p is the aver- 
age momentum). Now it is clear that the decrease of the 
temperature by the factor of 2 in experiment!^ might de- 
crease the normal density and therefore, shear viscosity 
at the center of the trap by a factor of 20. Although the 
effective ID viscosity should be obtained as a result of 
a complicated averaging of the two-fluid hydrodynamic 
equations in directions transverse to the trap this rough 
estimate shows that the possibility to observe the dis- 
persive behavior in quasi ID traps is indeed feasible and 
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requires the experiments at somewhat lower temperature 
than the one in Ref. EH 

To show that the presented analysis of the shock wave 
physics of KdVB equation describes qualitatively the so- 
lution of actual hydrodynamic system ( 4|5 1 we present 
the numerical solution of hydrodynamic equations in 
Fig. [3] for the hydrodynamic system derived for a uni- 
tary gas in a quasi-ld trap (see Ref. [TO] and the section 



VIB2|of this paper). 



We obtain numerical solutions of the continuity and 
Euler equation s using the technique of smoothed-particle 
hydrodynamic jl S I 11 ! 70 ! 71 !, This numerical method is a 
mesh-less algorithm. The algorithm essentially consists 
of mapping the system of differential equations describ- 
ing density and velocity fields to dynamic equations of a 
sufficiently large set of moving pseudo-particles (approx- 
imating the density and velocity fields) and then solving 
the molecular dynamics of these pseudo-particles. 

The results of a numeric solution of hydrodynamic 
equations are shown as a red (dashed) plot in Fig. 3] 
The parameters in equations are taken similar to the 
conditions of the experiment 10 discussed in the section 
|VI B 3| We used a constant background with the back- 
ground density given by Eq. |102| with the initial density 
profile of a Gaussian form of the height equal to po, ie, 
in Eq. |107| we take e = 1. The full- width at half max- 
ima of the Gaussian bump is taken to be about 20 fjxa 
(5% of system size) similar to the experiment)^. The red 
(dashed) plot in Fig. [3] corresponds to a v = 1 and we 



clearly see that the shock waves are of dissipative na- 
ture. Qualitatively the red (dashed) plot is very similar 
to the evolution described by the Burgers equation ( 35 1 . 



We also plot (black plot) the profile obtained numerically 
for a v = 0.1 (see the discussion of the feasibility of ob- 
taining small viscosities in finite systems above). For this 
value of effective dissipation we clearly see non- vanishing 
oscillations in the stationary profile of the shock wave. 



VII. CONCLUSION 

We presented a unified picture of dynamics of a large 
class of nonlinear systems in the limit of weak nonlin- 
earity, dissipation and dispersion (WNDD). This picture 
is delineated by the Korteweg-de Vries-Burgers equation 
( 25 1 for the deviation of the density from its background 
value. We described various regimes of nonlinear evolu- 



tion described by KdVB and summarized these regimes 
in the diagram in Figure [T] The parameters of KdVB 
have been related to those of a generic one-dimensional 
fluid and then to particular parameters of several models 
used in cold atom research. These relations are summa- 
rized in the Table [TJ 

Arguably the most interesting regime in WNDD limit 
is the formation of dispersive shock waves which are 
known to occur in interacting systems ^ 9 * 20 * 83 ! The anal- 
ysis presented in this paper uses the well-known results 
in the theory of nonlinear partial differential equations. 



It allows to make quick estimates of important scales and 
parameters of shock waves. 

Although in this work we focused on a Galilean invari- 
ant fluid, the general analysis can be extended to other 
systems as well. KdVB equation still remains the most 
universal description in WNDD limit unless additional 
symmetries restrict the form of the terms of the lowest 
order in nonlinearity and/or field gradients (e.g., the un- 
derlying particle-hole symmetry might require the sym- 
metry u — u so that the lowest allowed nonlinear term 
has a form of u 2 u x leading to the so-called modified KdV 
equation). The WNDD limit therefore allows to divide 
the complicated dynamics of nonlinear systems into sev- 
eral distinctive "nonlinear universality classes" giving an 
insight into their generic dynamic features (for a recent 
related discussion see Refs. IT51 1541 and IM]) . 

We have also discussed another important way of gen- 
eralizing KdVB description appropriate for systems with 
long-range interactions. For those systems the locality 
of the effective description is not required. The relax- 
ation of the locality requirement leads to a wider class 
of integro-differential effective equations (see Appendix 
[A] for some examples). 

Another important assumption used in this paper is 
the possibility to have a consistent effective descrip- 
tion of the system by a single-component fluid. Of 
course, the presence of several components such as nor- 
mal/superfluid components in BECs or equivalent com- 
ponents in spinor BECs can be crucial for the accurate 
description of corresponding physical systems. 

The described reduction of nonlinear three- 
dimensional dynamics to effective chiral one-dimensional 
equation obtained in WNDD limit (25 1 cannot capture 



all fascinating and complex effects in collective behavior 
of interacting 3d systems. For example, taking into 
account the next order terms in the reductive pertur- 
bation expansion will lead to the interaction between 
chiral sectors (left and right moving waves) neglected in 
this work. The next order terms in time derivatives will 
result in the frequency dependence of kinetic c oeffic ients 
essential for e.g., the physics of the unitary ga d 80 ! 85 !. 

Finally, the description used in this work is completely 
classical (for example, there is no Ti in Sees. II-V. The 
quantum physics enters here only through the values of 
effective parameters of classical hydrodynamic equations 
(see Sec. VI). One of the most interesting problem is to 
identify the quantum effects not describable by conven- 
tional classical physics (quantization of vorticity in su- 
perfluid system can serve as an example of such effects) . 
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Appendix A: Relaxing locality 

In the main part of the paper we assumed that the sys- 
tem of particles is well described by a local hydrodynam- 
ics theory. The assumption of locality was important in 
deriving the effective KdVB equation for those systems. 
In the presence of long range interactions between par- 
ticles our derivations should be modified and generally 
different counting scheme should be used in the reduc- 
tive perturbation method. In this appendix we show how 
the assumption of locality can be relaxed on the exam- 
ple of two models: the Calogero model and the dipolar 
BEC model. The former is often used as a theoretical toy 
model as it is exactly integrable and very well studied. 
The la tter is relevant to some experimental cold atom 
systems^EH. 



1. Calogero Model 

The Calogero model is the model of one-dimensional 
particles interacting via the inverse square potential 86 . 
The Calogero Hamiltonian is given by 



H = 



1 N 

— Y 



2m 

3 = 1 



2m 



N 

E 



A(A-l) 



Xk, 



(Al) 



Here the coupling constant A is dimensionless. For sim- 
plicity we consider here the quasiclassical limit corre- 
sponding to A 3> 1 and in the following replace A (A— 1) — > 
A 2 and neglect corrections of the order of 1/A. It turns 
out that one can write dow n the collective field theory of 
the above model exactly^ * 3 1 I 87 I, ie, including all gradi- 
ent and nonlinear terms. This field theory can be again 
casted in the form ( 4|5 ) with no dissipation r\B = and 



w(p) 



A{r) = 



n 2 \ 2 

m? 

h 2 x 2 

2m 2 



(npf 



^Px 



(A2) 
(A3) 



Here the superscript H stands for the Hilbert transform 
defined by 



f H (x)=*P.V.[ 



f(y) d v 



y 



(A4) 



Notice that the second term in ( A2 ) is nonlocal reflecting 



the long range interaction between Calogero particles in 

Eq.rni 



The sound velocity in the collective field theo ry for 
Calogero modeP^D can be obtained from (A2| using 



( 24 1 and is given by 



c = — Xtrpo ■ 
m 



(A5) 



The second term in (A2| is of the first order in d x and 
the Burgers counting scheme ( 36 39 1 should be used in 
the reductive perturbation method of Sec. |IV| As a result 
we obtain the so-called Benjamin-Ono equation 

u t + (uu£ + ctBO { u ii) H = (Benjamin-Ono) 

_(A6) 

as an effective description of Calogero hydrod ynam ics in 
the limit of weak nonlinearity and dispersion 30 ! 31 : The 
coefficients of ( A6 ) are given by 

(AT) 



rn 



Two remarks are in 



2m 

order. First. 



(A8) 

notice that the 



Benjamin-Ono term (the last term of Eq. A6l dominates 



the conventional KdV term in the long wave limit. 
It is of the same order as the Burgers term u^. Sec- 
ond, although the term scales similar to the Burg- 
ers term, it is very different in nature as it does not re- 
sult in any dissipation. Indeed, the Fourier transform 
(uf^)k = ik\k\uk and this term results in the real cor- 
rection lu ~ asok\k\ to the spectrum of linear waves in 
contrast to the Burgers term giving imaginary (dissipa- 
tive) correction lu ~ —ivk 2 . 



2. Quasi- ID Dipolar BEC 

Systems of bosons interacting via long range potential 
have been recently realized experimentally^ 3 -^ for bosons 
with dipole interactions. The quasi-lD dipolar BEC can 
be described by a non-local versioifHESES f ^ ne Q ross _ 
Pitaevskii 61 equation. By choosing a sufficiently large 
radial trap frequency it is possible to freeze the radial 
motioipSEQ of the BEC and thereby obtain a quasi- ID 
dipolar BEC described by the following non-local one- 
dimensional GPE (see Ref. [SO] for the derivation) 



iUdtip(x, t) = I — 



2m 



d 2 + mw( 



) UOM), (A9) 



where 



w{p) 



ml i 



(A 1DP + Bl 2 ± d 2 p u ^ 



(A10) 



Here g is the 3D contact interaction strength, lj_ z — 
\Jhj (mcu± iZ ) is the transverse (axial) oscillator length, 
and the dimensionless constants Aid and B are 



Aid = 
B = 



1 

2tt 
3 

8^ 



1 + \tdd (1 - in 2 ) 
£dd (l - in 2 z ) . 



(All) 
(A12) 
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n z is the z component of the direction of the dipole axis 
n and tdd is a natural dimensionless parameter for the 
relative strength of dipolar and s-wave interactional. 
In (A10) we used the dipolar integral transform 



fw = 

defined by its kernel 
U w (x) = 



U 1D {x - y)f(y)dy 



2l\ 



erfc 



y/2l ± 



(A13) 



(A14) 



where erfc is the complementary error function. The 
Fourier transform of U\d{x) is given by 



U 1D (k) 



dn - 



i K+(klj_y 
and its asymptotic behavior as k — ¥ is 

f(kl ± ) 



U 1D (k) ~ -j E - In 



(A15) 



(A16) 



where 7^ 0.5772 is the Euler-Mascheroni constant. 
The hydrodynamic c hang e of variables ip = 

•Jpe"* ->o v< - x ^ dx brings (A9) to the hydrodynamic 
forrrPS (4j5l with the chemical potential w(p) given by 
dAlOt and 



A{ P ) 



2m? 



(A17) 



The speed of sound in this system can be easily found 
using ( |24|A10| 



9Po 
ml 2 ^ 



A 



ID 



ml 



-(jA 



ID 



.1/2 



(A18) 



where we introduced the dimensionless coupling constant 



mpo 



(A19) 



The resulting Euler equation in this case is non-local. 
The perturbative scaling scheme (Eqns. 17|18|2l[22 l 
gives rise to the following non-local-KdV-like equation 
("dipolar-KdV") 



with 





(dipolar-KdV) 


(A20) 


-- 


(A21) 




h 2 




a = a 2 ' 


(A22) 


8m z c 


h 2 


(A23) 


a D =lB . 

2m z c 



Here the superscript U again denotes the transform 
( A13 1 . As we work only in the long wavelength limit we 



can think of the transform's kernel as of Fourier trans- 
form of (A16). Although, technically in the asymptotic 
long wavelength limit k — > the ac term will domi- 
nate the a term we keep both of them in (A2CM. In- 



deed, the kernel (Al6l is growing logarithmically slow 
with k and both terms become of the same order at 
kl± ~ exp(7ra/a£>) which is defined by the actual values 
of the dimensionless parameters of the dipolar system. 
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